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A cohesive element for shell analysis is presented. The element can be used to simulate 
the initiation and growth of delaminations between stacked, non-coincident layers of shell 
elements. The procedure to construct the element accounts for the thickness offset by 
applying the kinematic relations of shell deformation to transform the stiffness and internal 
force of a zero-thickness cohesive element such that interfacial continuity between the layers 
is enforced. The procedure is demonstrated by simulating the response and failure of the 
Mixed Mode Bending test and a skin-stiffener debond specimen. In addition, it is shown that 
stacks of shell elements can be used to create effective models to predict the inplane and 
delamination failure modes of thick components. The results indicate that simple shell 
models can retain many of the necessary predictive attributes of much more complex 3D 
models while providing the computational efficiency that is necessary for design. 
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coefficients of thermal expansion 
displacement jump for damage initiation 
final displacement jump 

= mode I, mode II, and mode III displacement jumps 
tangential displacement jump 

difference between the curing temperature and room temperature 
Kronecker delta 

equivalent displacement jump 
viscoelastic stabilization factor 

parameter of Benzeggagh-Kenane mixed-mode fracture criterion 

interfacial tractions 

nominal peel and shear strengths 
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I. Introduction 

Delamination is one of the predominant modes of failure in laminated composites when there is no 
reinforcement in the thickness direction. In order to design structures that are flaw- and damage-tolerant, analyses 
must be conducted to study the propagation of delaminations. The calculation of delaminations can be performed 
using cohesive elements 1 ' 2 , which combine aspects of strength-based analysis to predict the onset of damage at the 
interface, and fracture mechanics to predict the propagation of a delamination. 

Cohesive elements are normally designed to represent the separation at the zero-thickness interface between layers 
of three-dimensional elements (3D). Unfortunately, the softening laws that govern the debonding of the plies often 
result in extremely slow convergence rates for the solution procedure, and require very refined meshes in the region 
ahead of the crack tip. In addition, it is often difficult to model zero-thickness interfaces using conventional finite 
element pre-processors. 
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The typical requirement of extremely fine meshes ahead of the crack tip 1 3 ’ 4 has been recently relaxed by a 
procedure that artificially reduces the interfacial strength while keeping the fracture energy constant in order to 
increase the length of the cohesive zone. A longer cohesive zone allows the use of larger cohesive elements 3 . 
However, the difficulties associated with the generation of complex meshes composed of stacks of 3D elements 
joined together with zero-thickness cohesive elements often render the use of cohesive elements impractical outside 
research environments. 

Shell elements are more efficient for modeling thin-walled structures than 3D elements. Shell finite element 
analysis has been used in conjunction with the Virtual Crack Closure Technique (VCCT) by Wang et al. 5 6 to 
evaluate strain energy release rates for the damage tolerance analysis of skin-stiffener interfaces. Wang used a wall 
offset to move the nodes from the reference surfaces to a coincident location on the interface between the skin and 
the flange. Therefore, compatibility was ensured by equating the three translational degrees of freedom of both 
shells with constraint equations and by allowing the rotations to be free. These models were shown to require a 
fraction of the degrees of freedom than are needed for full 3D analyses while producing accurate values for mode I 
and mode II strain energy release rates. However, Glaessgen et al. 7 observed that shell VCCT models of debond 
specimens in which the adherends are of different thickness predict the correct total energy release rate but that the 
mode mixity does not converge with mesh refinement. This phenomenon is similar to the well known oscillatory 
behavior of the in-plane linear elastic stress field near the tip of a bi-material interfacial crack 8 . 

Very few references discuss the use of shell models in combination with cohesive elements and most of them are 
designed for explicit time integration. A shell model for discrete delaminations under mode I loading using finite- 
thickness volumetric elements connecting shell elements was proposed by Reedy 9 . Johnson et al. 10 applied the 
bilinear traction-displacement of the Crisfield 11 cohesive model in conjunction with the Ladeveze mesomodel 12 for 
intralaminar damage to investigate wing leading edge impact and crash response of composite structures. A mixed- 
mode adhesive contact shell element that accounts for the thickness offset and the rotational degrees of freedom in 
the shells was developed by Borg 13 for the explicit finite element code LS-DYNA®. Bruno et al. 14 developed shell 
models that use rigid links to offset the location of the cohesive element nodes from the shell reference surface to the 
interface between the layers. Using this technique, multi-layered shell models were analyzed using the zero- 
thickness cohesive elements available in the commercial finite element code LUSAS®. 
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Cohesive elements have been found useful to study fracture along bi-metallic interfaces 15 " 18 . Cohesive laws 
eliminate the stress singularity at the delamination front and, consequently, they prevent the oscillation in the stress 
field and the corresponding non-convergence of fracture mode mixity 18 . Therefore, it is expected that the non- 
convergence of the mode mixity observed by Glaessgen in shell VCCT models, which is also caused by a mismatch 
of the elastic stiffnesses of the adherends, is also avoided in cohesive models. However, a study of the mesh 
independence of mode mixity in cohesive shell models has not been attempted by the present authors. 

The objective of this work is to present a method for the construction of shell cohesive elements for implicit 
analysis. The formulation of an 8-node cohesive element for modeling delamination between shell elements on the 
basis of a 3D cohesive element previously developed by the authors 1, 2 is proposed. The formulation takes into 
account the rotational degrees of freedom at the shell nodes and the relative location of the interface which allows 
the construction of models with more than one layer of cohesive interfaces and does not require coincident nodes 
along an interface nor does it require the use of rigid links. Therefore, the formulation allows the use of several 
layers of shell/cohesive elements through the thickness of a thick laminate. Since the cohesive elements have the 
topology of a standard 8-node element, the models are simpler to generate than those based on zero-thickness 
interfaces. 

In the following sections, the constitutive damage model previously developed by the authors will be outlined. 
Then, a procedure to construct a shell cohesive element from that of any conventional zero-thickness cohesive 
element will be described. Finally, three different examples of increasing structural complexity will be presented to 
show that simple shell models based on the proposed cohesive element for shells can be used to represent the onset 
and propagation of delamination in structural components. 


II. Bilinear Cohesive Law 

A cohesive constitutive law relates the traction, r, to the displacement jumps, A, at the interface. The bilinear 
softening model, which is chosen here for its simplicity, is shown in Fig. 1. One characteristic of all softening 
models is that the cohesive zone can still transfer load after the onset of damage (A 0 in Fig. 1). For pure mode I, II or 
III loading, after the interfacial normal or shear tractions attain their respective interlaminar tensile or shear 
strengths, the stiffnesses are gradually reduced to zero. The area under the traction-displacement jump curves is the 
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respective (mode I, II or III) fracture energy. Using the definition of the J integral, it can be shown that for small 
cohesive zones 19 : 

| f r (A) dA. = G c (1) 

where G c is the total area under the traction-displacement law, is the critical energy release rate for a particular 
mode, and A-' is the displacement jump at failure, shown in Fig. 1. The penalty stiffness K is an arbitrarily large 
number selected such that the presence of undamaged cohesive elements does not introduce appreciable compliance 
to the structure 3 . 

A bilinear cohesive law is assumed for each of the three fracture modes. It is assumed that the 3 direction is 
normal to the interface and that the interlaminar shear strength T° shear is independent of the shearing direction. Then, 
the displacement jumps for damage initiation in each mode are simply 


Mode I: 


Mode II: 


A°3=^ 

3 K 


A 0 _ J shear 


K 


^0 ^ shear 

Mode III: ' K 

Similarly, the final displacement jumps are proportional to their corresponding toughness G ic as 
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A. Cohesive Law for Mixed-Mode Delamination 


In structural applications of composites, delamination growth is likely to occur under mixed-mode loading. 
Therefore, a general formulation for cohesive elements must deal with mixed-mode delamination growth problems. 

Under pure mode I, II or III loading, the onset of damage at the interface can be determined simply by 
comparing the tractions with their respective allowable. Under mixed-mode loading, however, damage onset may 
occur before any of the traction components involved reach their respective allowable. Therefore, a mixed-mode 
criterion must be established in terms of an interaction between components of the energy release rate. 

The constitutive damage model used here, formulated in the context of Damage Mechanics (DM), was 
previously proposed by the authors 2 . The damage model simulates delamination onset and delamination propagation 
using a single scalar variable, d, to track the damage at the interface under general loading conditions. An initiation 
criterion that evolves from the Benzeggagh-Kenane fracture criterion 20 (B-K) ensures that the model accounts for 
changes in the loading mode in a thermodynamically consistent way and avoids restoration of the cohesive state. 

The constitutive model prevents interpenetration of the faces of the crack during closing, and a Fracture 
Mechanics-based criterion is used to predict crack propagation. The parameter A is the norm of the displacement 
jump tensor (also called equivalent displacement jump norm), and it is used to compare different stages of the 
displacement jump state so that it is possible to define such concepts as ‘loading’, ‘unloading’ and ‘reloading’. The 
equivalent displacement jump is a non-negative and continuous function, defined as: 

A = ^(^3) +(\hear) 2 ( 4 ) 

where < • > is the MacAuley bracket defined as < x >= t(x + |x|), which sets any negative values to zero. The term A 3 

is the displacement jump in mode I, i.e., normal to the mid-plane, and A shear is the tangential displacement calculated 
as the Euclidean norm of the displacement jump in mode II and in mode III: 

A,w=a/(A,) 2 +(a 2 ) 2 (5) 

A bilinear cohesive law for mixed-mode delamination can be constructed by determining the initial damage 
threshold A n from the criterion for damage initiation and the final displacement jump, A ; , from the formulation of the 
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propagation surface or propagation criterion 2 . For the B-K fracture criterion, the mixed-mode displacement jump for 
damage initiation is 2 


A° = 




G n +G ni 

G T J 


(6) 


where the B-K parameter rj is obtained by curve-fitting the fracture toughness of mixed-mode tests and the fracture 
mode ratio is 2 
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and where the displacement jump ratio is defined as P = A sheal . /( A shear + (A 3 )). 

The displacement jump for final fracture is also obtained from the critical displacement jumps as 


A3A3 + (A sheal .A r shear A3A3 ) 


Pa/- 
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( 8 ) 


During overload, the state of damage d is a function of the current equivalent displacement jump A 


, a/(a-a°) 

The corresponding tractions can be written as 


( 9 ) 


r i = D ij^j = A'ijK 
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(10) 


where the Kronecker c) (/ is used to prevent the interpenetration of the surfaces of a damaged element when contact 
occurs 2 . 

In summary, by assuming that toughnesses and strengths in modes II and III are equal to each other, it is shown 
that the mixed-mode constitutive equations of a cohesive element are defined by five material properties and three 
displacement jumps: G Ic , G IIc , r® , T° shear , K, A b A 2 , and A 3 . This formulation of the damage model allows an explicit 
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integration of the constitutive model and ensures consistency in the evolution of damage during loading, unloading, 
and changes in mode mixity. The derivation and implementation of the cohesive damage model is described in Refs. 
[2, 21] and the reader is referred to those sources for details that will not be repeated here. 


B. Viscoelastic stabilization 

Material models exhibiting softening behavior and sharp stiffness degradation typically exhibit severe 
convergence difficulties in implicit analysis programs. To improve the convergence rate of the iterative procedure, a 
viscous stabilization scheme was implemented as suggested by ABAQUS 22 for its cohesive element and composite 
damage models. The procedure is a generalization of the Duvaut-Lions regularization model which, for sufficiently 
small time increments, maintains the tangent stiffness tensor of the softening material as positive definite. Defining 
d as the “kinematic” damage that depends on the displacement variables and d v as the “stabilized” damage variable, 
the rate of change in the damage variable d v is a function of the damping factor p: 


d r =— (d-d v ) 
P 

The stabilized damage variable at time t+At is updated as 


(11) 


- At (j 

1 P H 1 ' 

t+At P +At l+At 

i , u 

p+At 


( 12 ) 


It can be observed that the viscoelastic factor has the units of time and, therefore, its effect depends on the 
pseudo-time load proportionality factor which is selected arbitrarily by the user. In the following sections, the effect 
of the pseudo-time on the viscoelastic factor is eliminated by specifying it as a nondimensional value, p/t F where 
t F is the estimated pseudo-time at the maximum load. 

The response of the damaged material is evaluated using the stabilized damage variables. Using viscous 
stabilization with a viscosity parameter that is much smaller than the time increment can improve the rate of 
convergence in the presence of softening processes without compromising the accuracy of the results. However, the 
use of larger values of the stabilization parameter can cause a toughening of the material response and produce a 
stable propagation of delamination in problems with unstable delamination growth. 



C. Finite Element Implementation of Zero-Thickness 3D Element 


Zero-thickness cohesive elements as shown in Fig. 2a are used to simulate the response of the interface 
connecting two laminae of a composite laminate. 
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Fig. 2. a) 3D cohesive elements b) shell cohesive elements. 


The constitutive equation of zero-thickness cohesive elements is established in terms of displacement jumps and 
tractions across the interface. The stiffness matrix for conventional 3D cohesive elements is constructed using 
standard isoparametric linear Lagrangian interpolation functions. In a local coordinate frame with x-y-z aligned such 
that the z-direction corresponds to the shell normal, the displacement jumps between the top and the bottom faces of 
the element are 


^3 


W 


W 

A! 

> = • 

U 

. — 

U 

3^2. 

3 D 

V 

top 

V 


= B U 


3 D 


(13) 


bot 


where B is the matrix relating the element’s degrees of freedom U3D to the displacements between the top and 


bottom interfaces. The Principle of Virtual Work for 3D cohesive elements can be written as 2 


8V\ d JX((l-D)c)BdTU 3B =8Vl n \ySdA (14) 

where D is a diagonal matrix composed the damage state d at the interface. The term I represents the identity matrix, 
and C is the undamaged constitutive matrix whose diagonal is composed of the penalty stiffnesses K, and S is the 
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stress tensor. The stiffness matrix K 3/) of the element is easily identified in the left hand side of Eq. 14 as the 


integral over the area A of the element: 


k 3Z3 =JX((i-d)c)b^ ( 15 ) 

and the right hand side of Eq. 14 provides the element’s internal force: 

F 3fl =J 4 B r S«L4 (16) 

Because the stiffness K 3£) is constructed from the displacements U 3 /), only three degrees of freedom can be active 
at a node, i.e., u, v, and w. The integration of the terms in Eqs. 15 and 16 is performed numerically using either 
Newton-Cotes or Gaussian integration. The element proposed is implemented in ABAQUS 22 as a user-defined 
element. 


TIT. Formulation Shell Cohesive Element 

Shell cohesive elements have a finite non-zero thickness because they connect nodes on the reference surfaces of 
two shells, as shown in Fig. 2b. In addition to the u, v, and w translational degrees of freedom, the nodes of a shell 
element also possess three nodal rotations 6„ 6„ and 0-. For kinematic continuity before delamination propagation, 
all three translational degrees of freedom at the interface must be the same in both adjacent plies, as shown in Fig. 3. 


layer 

interface 


Fig. 3. Kinematic continuity is achieved when the displacements at the interface are equal in both plies 24 . 

When using shell elements based on First Order Shear Deformation, the three displacements at the interface are 
calculated from the nodal degrees of freedom as 23 : 
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A. Spatial Orientation of the Shell Normals 

The transformations shown in Eqs. 20 are only valid when the normals to the upper and the lower shells are 
aligned with the z-direction of the global reference frame. In the case of arbitrary orientation of the shells in space, it 
is necessary to define local reference frames at each corner of the shell cohesive elements. In the present model, the 
local normal direction (z') is defined by the vector between the lower and upper nodes at a corner. Using a rotation 
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tensor, the normal and tangential components of the displacement jump can be expressed in terms of the 
displacement field in global coordinates. 

B. Special Configurations 

It was shown in Eqs. 17 and 20 that the element formulation is based on the thickness t“ and t l of the upper and 
lower shells being connected together by a cohesive element. Therefore, the values of the thicknesses must be 
provided for each corner of every cohesive element in the model by using element property data cards. In models 
where the thicknesses of the adherends vary across the model, different property cards would have to be generated 
for every cohesive element. To improve the ease of generating a model where the shell thicknesses are non-uniform, 
it is convenient to define common geometric configurations that allow the use of single property cards for large 
areas of cohesive elements. Two such configurations are shown in Fig. 4. The configuration shown in Fig. 4a, which 
is typical of skin/ flange constructions, can be defined by specifying the thickness t L of the skin, which is constant for 
all cohesive elements in the feature. For configurations that are symmetric with respect to the interface, as in the 
example shown in Fig. 4b, the ratio t u lt l is constant. In either case, the two thickness t“ and t l can be determined 
within the element with the additional knowledge that the distance t“ + t l is equal to twice the distance between the 
upper and the lower nodes and a single property card can be used for all elements in the feature. 



a) constant 



Fig. 4. Special configurations with constant thickness properties. 
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IV. Example and Verification Problems 


A. Mixed Mode Bending (MMB) Specimen 

The MMB test offers the possibility of measuring different mixed-mode ratios, ranging from pure Modes I to II, 
by changing the length e of the loading lever shown in Fig. 5. The loading arm is modeled with a multi-point 
constraint equation between the displacements at the load application point and the two loading points on the upper 
arm of the specimen. 

The MMB test of a 24-ply unidirectional AS4/PEEK (APC2) carbon fiber reinforced composite was simulated 
with a simple model composed of two shell layers, and the results were compared with the results of a 3D model 
previously published by the authors 2 . The specimen is 102 mm-long, 25.4mm- wide, with two 1.56 mm-thick arms. 
The initial delamination length is 34mm, and the length e of the lever is 43 mm, which was calculated to provide 
a mode mixity G/Gf= 50%. The material properties of the specimen are shown in Table 1, and the properties of the 
interface are shown in Table 2. It should be noted that the parameter tj was calculated from a fit of experimental data 
corresponding to a single test specimen at each mode ratio 25 . 



Fig. 5. MMB test configuration and finite element models. 
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Table 1. Material Properties of AS4/Peek' \ 


En(GPa) E 22 -E 33 (GPa) Vi 2 =Vi 3 

V 23 Gl 2 = 

=Gi 3 (GPa) G 23 (GPa) 

122.7 10.1 0.25 

0.45 

5.5 3.7 

Table 2. Properties of the interface 1 ' 25 . 

G Ic (N/mm) G IIc (N/mm) 

r 3 ° (MPa) 

r lar (MPa) n 

0.969 1.719 

80 

100 2.28 


The experimental results shown in Fig. 6 were performed by Reeder 25 . They relate the load to the displacement 
of the point of application of the load F on the lever. The analytical solutions shown in the figure were calculated 
using the equations from Reeder 26 with G c =1.123 N/mm, which corresponds to the point where the delamination 
growth was first observed. This point is the second experimental data point in Fig. 6. Therefore, this point lies on the 
analytical constant-G curve. It can be observed that the test data indicates an R-curve effect consisting of an increase 
in toughness that ranges from an initiation value of 1.123 N/mm to a propagation value of approximately 1.27 
N/mm. 

The shell finite element model shown in Fig. 5 was developed using two layers of shell elements. Each layer is 
composed of a mesh of 2 X 150 elements. Therefore, the size of the elements in the direction of crack propagation is 
0.68mm. The shell model only represents 1/10 of the specimen width, so the predicted applied load was multiplied 
by a factor of 10. The loading lever was modeled as a rigid link using a multi-point equation. To improve the 

convergence rate, a viscoelastic factor of pi t F = 2.0 • 10 4 was used, where t F is the estimated pseudo-time at the 
maximum load. A penalty stiffness of K= 1 O' N/mm 3 was used. Recommendations on the determination of K are 
provided in Ref. 4. The corresponding load-displacement results are shown in Fig. 6. For comparison, the results of 
a simulation using a three-dimensional model by the authors 2 is also shown in Fig. 6. It can be observed that both 
models predict the analytical peak load and load-displacement curve accurately. The correlation of the analytical and 
finite element models with the test results could be improved using the toughness for propagation rather than the 
value for initiation. For an applied displacement larger than 7 mm, the delamination front proceeds under the mid- 
span load point and the present analytical solution is no longer valid. At that point, the finite element models diverge 
from the analytical solution. 

To verify the effects of mesh size and viscoelastic stabilization, a new shell model was created with 300 
elements in the propagation direction, or twice the number of elements of the base model. The model was analyzed 
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with p= 0. The predicted load with the refined mesh model is less than 1% lower than the load obtained from the 
previously discussed coarser model with plt F = 2. • 10 4 , as can be observed in Fig. 7. However, the refined mesh 
model failed to converge after a short amount of crack propagation. 

An alternative to using fine meshes of cohesive elements was proposed by Turon et al. ’ The approach consists of 
reducing the strengths and T° shear such that the length of the cohesive zone increases to cover at least three 
elements. According to Turon, the strengths that are required to achieve convergence using an element size of 0.68 
mm are: 


r 3 


J_ ( 9^33 G lc 
r?]l 32x5 l ele 


0.63 


^ shear . 1 1 9^^33 G 1Ic 

T ! Lar ^ shear V 32x5I ele 


0.67 


(21) 


Equations 21 indicate that the strength reduction needed for mesh convergence is greater for peel than for shear. 
Since both strengths should be reduced by the same factor to maintain a constant mode mixity, the reduced strengths 
r 3 and T shear that ensure that the cohesive zone is covered by at least five elements are: 
= 0.63 xr 3 ° = 50.4MPa ; r* /jear = 0.63 x r G shear = 63.0MPa . 



012345678 
Applied displacement, mm. 

Fig. 6. Predicted and experimental response of 50% mode mix MMB test. 
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Fig. 7. Effect of mesh density and viscoelastic stabilization on predicted load displacement response of 
50% mode mix MMB test. 

An analysis was performed with the original (coarse) mesh model with reduced strengths and no viscoelastic 
stabilization. The results shown in Fig. 7 indicate that the results from the original (coarser) mesh using reduced 
strengths are virtually identical to those using a refined mesh. In the present problem, the use of lowered strengths 
reduces the occurrence of divergence in the iterative solution and the use of viscoelastic stabilization is not required. 

B. Debonding of Skin-Stiffener Specimen 

Most composite components in aerospace structures are made of panels with co-cured or adhesively bonded 
frames and stiffeners. Testing of stiffened panels has shown that bond failure at the tip of the stiffener flange is an 
important and very likely failure mode. Comparatively simple specimens consisting of a stringer flange bonded onto 
a skin have been developed by Minguet et al." 7 to study skin/stiffener debonding. The configuration of the 
specimens shown in Fig. 8 was analyzed by Cvitkovich et al. 2S The specimens are 203 mm-long, 25.4 mm-wide. 
Both skin and flange were made from IM6/3501-6 graphite/ epoxy prepreg tape with a nominal ply thickness of 
0.188 mm. The skin lay-up consisting of 14 plies was [0/45/90/-45/45/-45/0]s, and the flange lay-up consisting of 10 
plies was [45/90/-45/0/90]s. 
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The properties of the unidirectional IM6/3501-6 unidirectional graphite/epoxy tape reported by Cvitkovich" 8 are 
shown in Table 3. However, the interface properties for IM6/3501-6 were not available. Instead, the interface 
properties of AS4/3501-6 are used, as suggested by Cvitkovich 28 . The values of G Ic , G IIc and the parameter 77=1.75 
for the B-K criterion shown in Table 4 are taken from test data for AS4/3501-6 29 . 


Tension Test 


203 mm 


Strain gages 


25.4 mm 


I I 


42 mm 
50 mm 


127 mm 





Ip 


I 1 

i 


Extensometer 

Fig. 8. Skin-stiffener specimen configuration 28 . 


Table 3. Material Properties of IM6/3501-6 unidirectional graphite/epoxy tape [Ref. 28] 


E„(GPa) E 22 =E 33 (GPa) 

Vi 2 — v [3 

v 23 

G 12 =Gj 3 (GPa) 

G 23 (GPa) 

144.7 9.65 

0.3 

0.45 

5.2 

3.4 

Table 4. Properties of the interface [Refs. 28,29,30]. 


G Ic (N/mm) G IIc (N/mm) 

r° (MPa) 

A°w (MPa) 

1 


0.0816 0.554 

61 

68 

1.75 



A complete analysis of the delamination growth in the tension specimen requires a high degree of complexity. 
For instance, Krueger 31 developed highly detailed two-dimensional models using up to four elements per ply 
thickness. In previously published work 30 , the authors determined that it is possible to predict the debond load of the 
specimen using cohesive elements using a relatively coarse three-dimensional model. To keep the modeling 
difficulties low and the approach applicable to larger problems, the 3D model was developed using only two brick 
elements through the thickness of the skin, and another two through the flange. The complete model consists of 
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1,002 three-dimensional 8-node C3DI elements and 15,212 degrees of freedom. This model did not contain any pre- 
existing delaminations. 

To account for the effect of residual thermal strains, a thermal analysis step was conducted before the 
mechanical loads were applied. The same coefficients of thermal expansion («;/=-2.4xlO" 8 /°C and a?2 = 3-7x 1 0" 5 /°C) 
are applied to the skin and the flange, and the temperature difference between the curing and room temperature is 
AT= -157 °C [Krueger 31 ]. The flange has more 90° plies than 0° plies and the skin is quasi-orthotropic, so that 
residual thermal stresses develop at room temperature at the skin/flange interface. 

In contrast to previous models of the skin-stiffener specimen, an extremely simple shell model was developed 
that uses one layer of shell elements to represent the skin and another layer to represent the flange of the stiffener. 
When subjected to tensile loads, the primary mechanism for debonding is shear lag. This fracture mechanism is 
captured in the model by using a layer of shell cohesive elements to bond the two shell layers together. The model is 
composed of six elements across the width, and the mesh in one of the two flange tip regions is refined. The element 
size in the propagation direction is 0.5 mm, and the rest of the flange is meshed with 1.0 mm elements. A penalty 
stiffness of K= 10 5 N/mm 3 was used for this analysis. The model is composed of 624 shell and cohesive elements 
and only 3192 degrees of freedom. Since this particular problem is dominated by the initiation of the debond rather 
than its propagation, the nominal strengths were not reduced, and an accurate calculation of the energy release rates 
during propagation was not attempted. 

A deformed plot of the finite element model immediately after flange separation is shown in Fig. 9. 



Fig. 9. Two-shell model of skin-stiffener specimen. 
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In previous work by the authors, it was stated that eigenmode analysis of the element stiffness matrices has 
shown that Gaussian integration can cause undesired spurious oscillations of the traction field when large traction 
gradients are present over an element and, consequently, that Newton-Cotes integration is superior to Gaussian 
integration in cohesive elements 1 ' 32 . However, in the analyses of the skin-stringer specimen, it was found that the 
Newton-Cotes integration causes odd wrinkling of the thin flange tip during the onset of delamination, as illustrated 
in Fig. 10a. The wrinkling deformations were found to delay the full separation of the flange tip. On the other hand, 
the wrinkling deformations were virtually eliminated by using Gaussian integration, as shown in Fig. 10b. 
Consequently, previously established conclusions on the preferred integration scheme may need to be re-examined. 



a) Newton-Cotes integration b) Gaussian integration 

Fig. 10. Onset of delamination with two different integration schemes (10X deformation). 


The predicted load-extensometer response for two values of viscoelastic stabilization are shown in Fig. 11. The 
initiation of delamination predicted using a viscoelastic stabilization factor of p/t F = 2.-1CT 4 is 17.2 kN, which is 
24% below the mean debond load measured by Cvitkovich 28 (22.7 kN). The unstable nature of delamination 
propagation can be deduced from the reduction of the applied load during the propagation of the delamination. The 
predicted delamination initiation load using p!t F = 2. - 10 3 is 17.5 kN. However, delamination propagation is 
stabilized by the use of a larger factor, as can be observed in Fig. 11. The advantage of using a higher amount of 
stabilization is that the solution is obtained much faster than with a lesser amount of stabilization: the model with a 
stabilization factor of 2. • 1 0 -4 took 6 minutes of CPU time on a 3.20GHz Xeon-based PC, while the model with a 
factor of 2. • 1 0 3 took less than 3 minutes. 
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The previously published results of a 3D model are also shown in Fig. 11. The slope of the force-extensometer 
curves for the 3D and shell models are different because nodal displacements were used to calculate the 
extensometer measurement. For the 3D model, the nodes are located on the surface of the laminates, while for the 
shell models, the nodes are located on the reference surfaces. It can be observed that the 3D model predicts a more 
accurate debond load of 21.9 kN, which is probably a result of the higher fidelity of the flange tip region in the 3D 
model. 



Fig. 11. Predicted and experimental load-extensometer response for skin-stiffener specimen in tension. 

C. Strength Analysis of Thick Composite Lug 

Shell finite elements are typically reserved for the analysis of thin-section structures such as fuselage skins. 
Since shell elements either neglect or simplify the transverse deformations, they are not usually applied to model 
components in which the thickness dimension is of the order of magnitude of other characteristic dimensions of the 
component. Flowever, computationally efficient models of thick components have occasionally been developed by 
stacking layers of shell elements. For instance, the complex through-the-thickness deformations of thick sections of 
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the Composite Armored Vehicle, which includes layers of glass/epoxy, ceramic tiles, and rubber, were successfully 
modeled by tying four or more layers of shell elements using multi-point constraint equations 24 . 

In addition to the computational efficiency that results from using shell elements, the use of shell cohesive 
elements to model layered structures provides the added benefit of simplifying the construction of complex models 
in several ways. Firstly, shell cohesive elements eliminate the laborious task of constructing multi-point constraints 
to tie the layers together. Secondly, unlike conventional cohesive elements, shell cohesive elements have a non-zero 
thickness and they have the same topology and connectivity as a standard volume element. Therefore, they can be 
constructed using any finite element pre-processor. Finally, shell cohesion elements do not require the use of wall 
offsets. Therefore, the construction of shell models is simplified and, more importantly, the shell reference surfaces 
do not have to be shifted to a common interface, which allows stacks of multiple shell layers. 

During the course of a failure analysis of a composite aircraft fin undertaken at NASA Langley 33 , a layered shell 
element of the critical attachment lug area was developed. The model detail shown in Fig. 12 illustrates how 
cohesive elements for shells can be used to construct models of large and geometrically complex thick composite 
structures by stacking layers of shell elements. The thick lug fittings were modeled with 14 layers of shell elements, 
which were connected with shell cohesive elements. All other regions of the model were modeled with a single layer 
of shell elements. The pin assembly was modeled as a rigid surface with a diameter equal to that of the lug hole. 
Frictionless contact equations were prescribed between the edge of the shells around the bolt hole and the rigid 
surface. This model has 20886 nodes with 34524 elements. The response of the layered-shell model was compared 
to a 3D model (no delamination) and it was found that the strain distributions predicted by the two models were in 
close agreement 33 . 
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Analysis with progressive 
delamination of inner fitting. 


Fig. 12. Layered shell model of composite lug with ply damage and cohesive elements. 

In addition to delamination at any interface between the 14 layers of shells, the model uses a progressive damage 
model in which the matrix and fiber damage is simulated by degrading the material properties. An example of 
damage analysis is shown in Fig. 12, where the regions in red represent areas where damage is predicted. The 
analyses for intra-ply damage were developed within ABAQUS using a USDFLD user-written subroutine. The finite 
element program calls this subroutine at all integration points of elements that have material properties defined in 
terms of the field variables. The subroutine call provides access points to a number of variables such as stresses, 
strains, material orientation, current load step, and material name, all of which can be used to compute the field 
variables. Stresses and strains are calculated at each incremental load step, and evaluated by the failure criteria to 
determine the occurrence of failure and the mode of failure. 

The results of a failure analysis indicate that a cleavage-type failure of the lug, in which an initial fracture 
develops in a direction parallel to the orientation with the greatest number of fibers. This direction is nearly aligned 
with the load direction, as shown in Figs. 13. 
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Fig. 13. Progressive failure analysis predicts a two-step cleavage failure of the composite lug 33 . 

The greatest advantage of the layered shell model is its ability to predict the initiation of delamination, as well as 
the component’s tolerance to delamination flaws. A study was undertaken to evaluate the damage tolerance of the 
lug in the presence of a pre-existing delamination. The flaw was introduced in the regions of the model shown in red 
in Figs. 14 by setting the region of shell cohesive elements shown in red to “failed” by defining the damage variable 
of these elements as unity. Failed cohesive elements do not have stiffness in tension or shear, but they prevent 
interpenetration of the crack surfaces in compression. 



Fig. 14. Configuration of delaminated lug and model for damage tolerance analysis. 

The strength analyses with and without a pre-existing delamination were performed under tensile loading. For 
this loading condition, the analyses did not predict delamination initiation nor propagation of the delaminated 
region. The load-displacement results of the two analyses are shown in Fig. 15. It can be observed that the initial 
flaw does not appreciably affect the overall stiffness or strength of the composite lug. The predicted loads for both 
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simulations fall between the load measured for the subcomponent test performed in 2003 and the W375 tension load 
that is believed to have caused the in-flight failure of the lug 33 . 



Fig. 15. Predicted lug strength with and without initial flaw. 


V. Conclusions 

The formulation of a cohesive element for shell analysis that can be used to simulate the initiation and growth of 
delaminations between stacked, non-coincident layers of shell elements was presented. The procedure to construct 
the element accounts for the thickness offset by applying the kinematic relations of shell deformation to transform 
the stiffness and internal force of a zero-thickness cohesive element such that interfacial continuity between shell 
layers is enforced. The internal load vector and the stiffness matrix of a cohesive element connecting shells can be 
calculated in a very simple way by means of a transformation matrix, independently of the cohesive constitutive 
models. 

Stacks of shell elements can be used to create effective models that predict the inplane and delamination failure 
modes of thick components. The formulation allows the use of multiple layers of shells through the thickness and it 
does not require offsetting nodes to a common reference surface. In addition, shell cohesive elements have the same 
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topology as three-dimensional elements so they can be generated using any finite element pre -processor while zero- 


thickness cohesive elements for three-dimensional elements cannot be generated with most pre-processors. 

Using representative examples of increasing structural complexity, it was shown that the proposed cohesive 
element for shells can be used to represent the propagation of a pre-existing delamination, as well as the onset and 
propagation of delamination in components that do not contain pre-existing cracks. The analysis of a 50% mode-mix 
MMB test indicates that the initiation and propagation of delamination under mixed mode can be predicted 
accurately. It was found that the convergence difficulties that are typical of cohesive laws can be mitigated with the 
use of viscoelastic stabilization or by reducing the interfacial strength to enlarge the process zone. The analysis 
results of a two-layer model of a skin-stiffener debond indicate that such a simplified model can produce results that 
may be sufficiently accurate for design and down-select analyses. For better predictive accuracy, models in which 
the flange taper is represented more accurately would be required. Finally, the analysis results of a large fin lug 
component were presented to illustrate how layered shell models with cohesive elements can be used to predict the 
strength of thick composite structures. This combination of cohesive elements and shells represents an effective way 
of simulating complex aeronautical structures that incorporate other non-linearities such as ply damage propagation, 
contact, and finite rotations. However, for the loading conditions that were analyzed, the cohesive elements in the 
model were not subjected to loads sufficiently large to cause any delamination growth. 
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